########
# What does this do?
#######
rm(list=ls())

# load libraries
packs = c('tidyverse', 'effects', "stargazer",
          'hrbrthemes', 'mlogit', "here")
source("r/LoadPkg.R")
loadPkg(packs)

# read data
covars = read_csv(here('repFile', 'data', 'covars_final.csv'))

# create outcome variables
## 3-level dv
covars$vote = recode(covars$vote_choice, 
                     `En contra` = 0L, 
                     `Did not vote` = 1L,
                     `A favor` = 2L)


# dummy uribe
covars$uribe = ifelse(covars$uribe_feel_fix == 'Favorable', 1, 0)
covars$santos = ifelse(covars$santos_feel_fix == 'Favorable', 1, 0)


# modeling
modDat =
  covars %>%
  dplyr::select(uribe, nse, age_factor, vote,
                voted_5_years_num, victim_num,
                woman, rearing_agg, cod_dane)


# modeling with mlogit()
mdata = mlogit.data(modDat[,c('vote', 'rearing_agg', 'nse',
                              'age_factor', 'voted_5_years_num',
                              'victim_num', 'woman', 'uribe')], shape = 'wide',
                    choice = 'vote')

m1 = mlogit(vote ~ 0 | rearing_agg + nse + age_factor + woman + uribe +
              victim_num + voted_5_years_num, data=mdata, reflevel = '2')

m2 = mlogit(vote ~ 0 | rearing_agg + nse + age_factor + woman + uribe +
              victim_num + voted_5_years_num, data=mdata, reflevel = '2', 
            alt.subset = c('2', '0'))

m3 = mlogit(vote ~ 0 | rearing_agg + nse + age_factor + woman + uribe +
              victim_num, data=mdata[mdata$voted_5_years_num==1,], reflevel = '2')

m4 = mlogit(vote ~ 0 | rearing_agg + nse + age_factor + woman +
              victim_num + voted_5_years_num, data=mdata, reflevel = '2')

# IIA test
hmftest(m1, m2)

# # clean up coef names
clean = names(coef(m1))
clean = str_replace_all(clean, pattern = '0:', replacement = 'Vote no: ')
clean = str_replace_all(clean, pattern = '1:', replacement = 'Abstain: ')
clean = str_replace_all(clean, pattern = 'rearing_agg', replacement = 'child-rearing')
clean = str_replace_all(clean, pattern = 'nseEstrato', replacement = 'NSE')
clean = str_replace_all(clean, pattern = 'nseEstrato', replacement = 'NSE')
clean[clean == 'Vote no: nseNo aplica (Sin estratificación)'] <- "Vote no: NSE = NA"
clean[clean == 'Abstain: nseNo aplica (Sin estratificación)'] <- "Abstain: NSE = NA"
clean = str_replace_all(clean, pattern = 'De',replacement = '')
clean = str_replace_all(clean, pattern = 'años', replacement = '')
clean = str_replace_all(clean, pattern = 'Más de', replacement = '')
clean = str_replace_all(clean, pattern = " a ", replacement = '-')
clean = str_replace_all(clean, pattern = "age_factor", replacement = '')
clean = str_replace_all(clean, pattern = "woman", replacement = 'Woman')
clean = str_replace_all(clean, pattern = "uribe", replacement = 'Uribe')
clean = str_replace_all(clean, pattern = "victim_num", replacement = 'Victim')
clean = str_replace_all(clean, pattern = "voted_5_years_num", 
                        replacement = 'Recent voter')


# table output
stargazer(m1, m3, m4, 
          covariate.labels = clean,
          keep.stat = 'n', style = 'io',
          keep = c(":rearing_agg", "intercept"),
          dep.var.labels = 'Voting status',
          title = 'Results from multinomial logit model. Reference category is having voted Yes referendum. Controls: Income, age, sex, Uribe supporter, victim, voted in a prior election.',
          label = 'mlogit',
          column.labels = c('Full sample', 'Recent voters', 'No Uribe control'),
          out = here("repFile", "paper", "figures", "mlogit.tex"))


stargazer(m1, m3, m4, 
          covariate.labels = clean,
          keep.stat = 'n', style = 'io',
          keep = c(":rearing_agg", "intercept"),
          dep.var.labels = 'Voting status',
          title = 'Results from multinomial logit model. Reference category is having voted Yes referendum. Controls: Income, age, sex, Uribe supporter, victim, voted in a prior election.',
          label = 'mlogit',
          column.labels = c('Full sample', 'Recent voters', 'No Uribe control'),
          out = here("repFile", "JOP submission", "mlogit.tex"))



#### multinomial effects
## convert factors to numeric for ease
mdata = 
  covars %>% 
  filter(nse != 'No aplica (Sin estratificación)')

# Convert to numeric
mdata$nse = 
  mdata$nse %>% 
  as_factor() %>% 
  fct_relevel("Estrato 1", "Estrato 2", 
              "Estrato 3", "Estrato 4",
              "Estrato 5", "Estrato 6") %>% 
  as.numeric()

mdata$age_factor = 
  mdata$age_factor %>% 
  as_factor() %>% 
  fct_relevel("De 18 a 25 años", 
              "De 26 a 35 años",
              "De 36 a 45 años",
              "De 46 a 55 años",
              "De 56 a 65 años",
              "Más de 65 años") %>% 
  as.numeric()


# Set parameters for model
dv='vote'
ivs=c('rearing_agg', 'nse', 'age_factor', 'woman', 
      'nse', 'uribe', 'victim_num',
      'voted_5_years_num')
form = formula(paste0(
  dv, ' ~ ', paste0(ivs, collapse=' + ')
) )

# relevel DV so baseline = 2
mdata$vote = 
  mdata$vote %>% 
  as.character() %>% 
  as_factor %>% 
  relevel(mdata$vote, ref = '2')

# Run model
# DV
multMod = nnet::multinom(form, data=mdata, Hess = TRUE)

# Get relevant parameter estimates
coefs = t( summary(multMod)$coefficients )
serrors = t( summary(multMod)$standard.errors )
zscores = coefs/serrors
pvals = (1 - pnorm(abs(zscores), 0, 1)) * 2

# get coefficient estimates
beta = multMod$wts[multMod$wts!=0]

# Set up scenario
rearRange = seq(min(mdata$rearing_agg), max(mdata$rearing_agg))
cov =
  mdata %>%
  dplyr::select(one_of(c('nse', 'age_factor', 'woman', 
                  'nse', 'uribe', 'victim_num',
                  'voted_5_years_num'))) %>%
  mutate(nse = as.numeric(as.factor(nse)),
         age_factor = as.numeric(as.factor(age_factor))) %>%
  summarise_all(funs(mean(., na.rm = TRUE)))

fake = cbind(1, rearRange, cov)

# Pull some draws
sims=1000
set.seed(6886); draws = MASS::mvrnorm(sims, beta, solve(multMod$Hess) )

# Calculate parameters
noPred = exp( draws[ ,str_detect(colnames(draws), '0:')] %*% t(fake) )
absPred = exp( draws[ ,str_detect(colnames(draws), '1:')] %*% t(fake) )

# Get denominator
simDenom = (1 + noPred + absPred )

# Get simulated probabilities for each category
genProb = 1/simDenom
noProb = noPred / simDenom
absProb = absPred / simDenom

# Pull out mean and 95% interval
info = function(x){ c( mean(x), quantile(x, probs=c(0.025, 0.975)) )  }
genProbSumm = t( apply(genProb, 2, info) )
noProbSumm = t( apply(noProb, 2, info) )
absProbSumm = t( apply(absProb, 2, info) )

# Set up dataframe for plotting
ggData = data.frame( rbind(
  cbind(rearRange, genProbSumm),
  cbind(rearRange, noProbSumm),
  cbind(rearRange, absProbSumm)
) )
colnames(ggData) = c('rearRange', 'mu', 'lo', 'hi')
ggData$scen = rep( c('Vote Yes', 'Vote No', 'Abstain'), each=nrow(fake) )

# Plot
ggplot(ggData, aes(x=rearRange, y=mu, ymin=lo, ymax=hi,
                       color=scen, fill=scen, shape = scen)) + 
  geom_ribbon(alpha=.2) + geom_line(size = 1) + 
  geom_point(size = 3) + 
  theme_ipsum_rc(grid = 'Y') + 
  theme(legend.position = 'top') + 
  labs(x = 'Child-rearing scale', 
       y = 'Predicted probability of outcome', 
       color = 'Vote status:', 
       fill = 'Vote status:', 
       shape = "Vote status:") + 
  scale_fill_manual(values = wesanderson::wes_palette(name = 'Darjeeling1', 
                                                      n = 3)) + 
  scale_color_manual(values = wesanderson::wes_palette(name = 'Darjeeling1', 
                                                       n = 3)) +
  scale_y_percent(limits = c(0, .7), breaks = seq(0, 1, by = .1))

ggsave(here("repFile", "paper", "figures", "multinom-probs.pdf"), 
       device = cairo_pdf)

ggsave(here("repFile", "JOP submission", "multinom-probs.eps"), 
       device = cairo_ps)

ggsave(here("repFile", "JOP submission", "multinom-probs.pdf"), 
       device = cairo_pdf)
